## 

rm(list = ls())

library(tidyverse)
library(pbapply)
library(rdrobust)
library(fixest)
library(broom)

## Def covars

covs <- c('turnout_party_2009_btw', 'soz_vers_beschaeftigte_share',
          'pop_density_km2', 'migration_out_share')

## Get federal election data

bt <- readRDS('data/data_federal.rds') %>%
  filter(!is.na(treated)) 

## List of outcomes

outcomes <- c("turnout_party", 
              "agg_left_party", 
              "agg_center_party",
              "agg_right_party")

## Source helper function to tidy rdrobust output

source("code/tidy_rd.R")

## Drop the 09 obs

bt <- bt %>% 
  filter(year > 2012)

## First differences

diff_df <- pblapply(outcomes, function(o) {
  out <- bt %>%
    filter(year > 2012) %>%
    pivot_wider(values_from = o, names_from = 'year', id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2017  - o2013) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(bt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs), state_id) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 0)

## Estimate for all outcomes

## Optimal BWs estimated separately

out_rd_opt <- pblapply(outcomes, function(o) {
  
  out <- rdrobust(y = diff_df %>% pull(!!o), 
                  x = diff_df$runvar, c = 0,
                  covs = diff_df[, covs])
  
  
  ## Return this : Robust B-C SE
  out2 <- out %>% 
    tidy_rd(se_nr = 3) %>% 
    mutate(outcome = o, method = 'robust bias-corrected') %>% 
    dplyr::rename(se = std.error, pval = p.value, 
                  bw = bw_left_h,
                  bw_bias = bw_left_b) %>% 
    mutate(conf.low90 = estimate - qnorm(0.95)*se,
           conf.high90 = estimate + qnorm(0.95)*se)
  
  ## Both
  
  out2
}) %>%
  reduce(rbind)

## Check

out_rd_opt %>%
  dplyr::select(outcome, method, estimate, pval)

##

out_rd_opt <- out_rd_opt %>%
  mutate(outcome = dplyr::recode(outcome, 
                                 `agg_left_party` = 'Left-wing\nparties',
                                 `agg_center_party` = 'Centrist\nparties',
                                 `agg_right_party` = 'Right-wing\nparties',
                                 `turnout_party` = 'Turnout')) %>% 
  mutate(outcome = factor(outcome, levels = unique(outcome)[c(2:4, 1)]))

# Figure A.13: Effect in R-P ----

p1 <- ggplot(out_rd_opt, 
             aes(outcome, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted',
             color = 'grey40')  +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.65) +
  geom_errorbar(aes(ymin = conf.low90, 
                    ymax = conf.high90), 
                width = 0, size = 1.05) +
  geom_point(shape =21,
             fill = 'white', size = 2.5) +
  xlab('') + ylab('RD estimate (p.p.)') + theme_bw()
p1
